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We review progress in the hydrodynamic description of heavy-ion collisions, focusing on 
recent developments in modeling the fluctuating initial state and event-by-event viscous 
hydrodynamic simulations. We discuss how hydrodynamics can be used to extract in- 
formation on fundamental properties of quantum-chromo-dynamics from experimental 
data, and review successes and challenges of the hydrodynamic framework. 



> 

1. Introduction 

G\ 

00 The large elliptic flow measured at the Relativistic Heavy-Ion Collider (RHIC) 

at Brookhaven National Laboratory and recently at the Large Hadron Collider 

t-H (LHC) at CERN is one of the most striking observations in heavy-ion collision 

experiments. This asymmetry of particle production in the transverse plane of the 
collision is interpreted as a sign of the system's hydrodynamic response to the 
initial geometry. Ideal hydrodynamics, which includes no viscous effects, was first 

. — used to describe heavy-ion collisions, and did a surprisingly good job in reproducing 

experimental data. Improved agreement with the data can be achieved using viscous 
hydrodynamic simulations with a very low shear viscosity to entropy density ratio. 
The applicability of hydrodynamics demands a short mean free path with respect 
to the system size. Therefore it is concluded that the created quark-gluon plasma 
is strongly interacting and behaves like a nearly perfect fluid. 

Recently interest has cascaded to higher harmonic flow, representing e.g. a tri- 
angular variation in the produced particle spectrum, which is non-zero in single 
events with lumpy initial energy density distributions. Hydrodynamic simulations 
together with models for the fluctuating initial state have been very successful in 
reproducing experimental data on all measured flow harmonics, their probability 
distributions and other quantities. 

The aim in the application of hydrodynamics to heavy-ion collisions and compar- 
ison to experimental data is the extraction of properties of the created quark-gluon 
plasma and to learn about the initial state and its fluctuations. This provides infor- 
mation on fundamental properties of quantum chromo dynamic systems. Here, we 
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review the tremendous progress that has been made over the last couple of years to 
achieve these goals. 



2. A short history 

We begin by reviewing the history of flow measurements at RHIC and LHC and the 
evolution of theory, in particular the hydrodynamic description of the bulk medium 
alongside experimental advances. 

It was widely expected before the first RHIC results in 2000 that, because of the 
extremely high energy densities achieved in heavy-ion collisions at ^/s = 200 GeV 
and the asymptotic freedom of QCD, the system would be one of weakly interacting 
partons. This system should behave like a gas and expand isotropically. However, the 
first results from RHIC demonstrated that this was not the case: Particles emerging 
from heavy-ion collisions at RHIC showed a significant azimuthal anisotropy. This 
anisotropy, characterized by the second Fourier coefficient w 2 in an expansion of 
the azimuthal particle distribution, was 50% larger than that measured at SPS 
(y/s = 17.3 GeV), and agreed well with calculations using ideal hydrodynamics, 
initialized to reproduce the measured multiplicity and transverse momentum spectra 
[I|j3j. The success of ideal hydrodynamics implied that the system was strongly 
interacting and very different from early expectations. 

Both measurements and theoretical calculations became more sophisticated. De- 
tailed measurements of vi for different particle species followed, showing a clear 
mass-splitting, predicted by hydrodynamics and due to the emergence of all parti- 
cles from a single velocity field. 

The relevance of viscous corrections for the hydrodynamic description of heavy- 
ion collisions was never completely ignored, but it was not until the ratio of viscosity 
over entropy density in strongly-coupled systems was calculated using AdS/CFT 
techniques 4 in 2004, in a regime where standard kinetic theory was known to break 
down, that the study of viscous effects became an urgent matter. The AdS/CFT 
calculations predicted that the viscosity to entropy density ratio of an TV = 4 super 
Yang-Mills quantum system in the strong coupling limit is rj/s = 1/47T, a value 
that can be rationalized by the argument that excitations cannot be localized with 
a precision smaller than their thermal wavelength [5j, but which had not been 
reliably calculated previously. 

From then on viscous hydrodynamic calculations and comparison to experimen- 
tal data emerged as the potentially best way to determine the transport properties 
of the matter created in heavy-ion collisions. Earlier perturbative QCD calculations 
predicted a rather large value of rj/s (compared to l/47r) with, however, significant 
error bands, mainly resulting from uncertainties in the relevant scales. [6] Non- 
perturbative lattice QCD calculations existed [7jj9], but because transport proper- 
ties are dynamical quantities, some doubts remain concerning their reliability. The 



quest to determine the effective rj/s with precision began and recent studies 10 



find it in the range 0.07 < rj/s < 0.43. This result is for an average value that corre- 
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sponds to the viscosity at the typical time when flow is built up. lOfll Some effort 



has been directed at developing a framework for precision determination of different 
parameters, like this effective shear viscosity to entropy density ratio, in a detailed 



multi-parameter \ 2 fit to multiple experimental quantities 12 . There are, however, 
missing ingredients like the fluctuations described below such that it is probably 
too early to attempt precision fits to experimental data. Future goals along these 
lines naturally include the determination of the temperature dependence of ry/s and 
other parameters, such as bulk viscosity. We will discuss these endeavors in more 
detail in the following sections. 

While theoretical work was focused on the study of viscous effects, experiments 
discovered the importance of fluctuations in the initial geometric configuration while 
studying Cu+Cu collisions. Copper ions were initially used in 2005 to study whether 
interesting physics effects found in Au+Au collisions would turn off. Interestingly, 



V2 was found to be quite large, even in central Cu+Cu collisions 16 , where i»2 
should be very close to zero if the initial energy density could be described as a 
convolution of two smooth nucleon distributions at impact parameter b « Ofm. 
The concept of the "participant eccentricity" , where the shape of the overlap region 
was not calculated relative to the classical impact parameter, but relative to an 
axis determined by the fluctuating participants, allowed for an explanation of the 
observed large i^- 

Fluctuations are also the reason why odd flow harmonics (such as V3, 1*5, etc.) are 
not zero. It was not until 2010 that real interest in the study of higher flow harmonics 



was sparked 17 18 , even though their potential was pointed out earlier [19]. First 
results on higher-order harmonic flow by the RHIC and LHC experiments appeared 
just before and at the Quark Matter 2011 conference. v\ through vq were shown to be 
sizable, each having its own amplitude and event plane angle. Odd flow harmonics 
were found to have a weak centrality dependence, characteristic of initial state 
geometric fluctuations. Furthermore, the v n become weaker with increasing n, as 
expected from the presence of viscosity, which more efficiently damps out higher 
order (smaller wavelength) fluctuations. 

The current state-of-the-art viscous hydrodynamic calculations with the most 
advanced models for the fluctuating initial state show remarkable agreement with 



flow measurements from both RHIC and LHC 20 . Even the event-by-event distri- 
butions of V2-V4 as measured by the ATLAS collaboration 21 are reproduced by 
these calculations. 

The historical developments outlined in this section are summarized in a time 
line in Fig. [I] presenting key figures and the increasing precision in the determina- 
tion of an effective shear viscosity to entropy density ratio of hot and dense nuclear 
matter and an outlook to future goals, including the determination of the tem- 
perature dependence of shear and bulk viscosity as well as more details on other 
transport parameters. The center column lists important developments as a function 
of time, including an outlook towards potential achievements in the next decade. 
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Important experimental and theoretical developments 
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Fig. 1. Time line of important experimental and theoretical developments leading towards in- 
creasingly precise understanding of flow, transport properties of the quark-gluon plasma, and the 
initial state and its fluctuations. On the left are three key figures [13})15] depicting the progress of 
hydrodynamic calculations and their success in describing experimental data, followed by a sketch 
of the uncertainty in a temperature dependent (jj/s)(T). On the right, the increasing precision 
in one key observable, the shear viscosity to entropy density ratio r)/s near its minimal value, 
is illustrated (see text for details). Figure adapted from the "Hot & Dense QCD White Paper", 
solicited by the NSAC subcommittee on Nuclear Physics funding in the US. 
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In the left column, we first highlight the early success of hydrodynamics in describ- 



ing results on transverse momentum dependent elliptic flow from RHIC 13 . The 
simulation described the data very well despite the fact that important ingredients 
such as a QCD equation of state and viscous effects were neglected. This early suc- 
cess was however reason to pursue further developments within the hydrodynamic 
framework, leading to more and more sophisticated modeling. One such important 
step is represented by the second figure, showing the comparison of first viscous 
hydrodynamic calculations with different initial conditions and different shear vis- 
cosity to entropy density ratios to experimental data. The third figure shows the 
comparison of experimental data on V2(j>t), "^(pt), ^4(pt), and v§(pt) from the 



ATLAS collaboration 22 to results from an event-by-event viscous hydrodynamic 



simulation including a QCD based model for the initial state fluctuations and us- 



ing a constant r)/s — 0.2, showing striking agreement 20 . The final figure is a 
sketch of the uncertainty in the temperature dependent value of T)/s, indicating 
that a combined analysis of LHC data and RHIC data at different energies using 
state-of-the-art hydrodynamic simulations could help determine this temperature 
dependence as well as that of other transport parameters. On the right, the increas- 
ing precision in one key observable, the shear viscosity to entropy density ratio rj/s 
near its minimal value, is illustrated. Shown results are from perturbative QCD 
calculations at leading logarithmic order 1 23 26 (shown formula) , complete leading 
order [6] (band from varying the scale by 20%), Anti-deSitter gravity/Conformal 
Field Theory (AdS/CFT) correspondence [3], lattice QCD - pure glue at ~ 1.6 T c , 



1.24 T c , and 1.58 T c , respectively [7H9 , ideal hydrodynamics 27 28 , perturbative 



QCD/kinctic theory 29 30 , and viscous hydrodynamics constrained by flow mea- 
surements [l4j 31-33 . We indicate that going forward, efforts should focus on the 
determination of temperature dependent quantities since more precise values of an 
effective rj/s depend on the collision energy and other parameters of the collision. 

Despite the great success of relativistic viscous hydrodynamics in describing 
experimental data from heavy-ion experiments, several puzzles remain. It is for 
example still not understood how the system can thermalize quickly enough to 
allow for the early applicability of hydrodynamics (which is needed to describe the 
data), or whether this early thermalization is even needed at all. We will discuss 
this and other open questions and issues in the remainder of this review. 



3. Theoretical framework 

The current standard for the viscous hydrodynamic description of relativistic heavy- 
ion collisions has been established with the derivation of hydrodynamic equations 
including all terms up to second order in gradients for a conformal fluid [34] . Addi- 
tional terms for non-conformal fluids with non-zero bulk viscosity have been derived 
subsequently [35] . 

In the ideal (non- viscous) case, the evolution of the system created in relativistic 
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heavy-ion collisions is described by the following 5 conservation equations 

W = o, « = o, (i) 

where T.^ is the energy-momentum tensor and Jg is the net baryon current. These 
are usually re-expressed using the time-like flow four-vector u M as 

Tff = (e + V)uW-Vgi", J£ = PB u», (2) 

where e is the energy density, V is the pressure, is the baryon density and 
g^ v = diag(l, — 1, — 1, — 1) is the metric tensor. The equations are then closed by 
adding the equilibrium equation of state 

V = V(e, PB ) (3) 

as a local constraint on the variables. 

In the first-order, or Navier-Stokes formalism for viscous hydrodynamics, the 
stress-energy tensor is decomposed into T^ t = Tff + , where Tff is given in 
Eq. ^ and the viscous part of the stress energy tensor is 

= r) (v"u v + VV - ^""Vo^j , (4) 

where = - u^u v is the local 3-metric and V^ 1 = A^<9„ is the local spa- 
tial derivative. Note that S^" is transverse with respect to the flow velocity since 
A M1/ u„ = and u v u v = 1. Hence, is a (timelike) eigenvector of the whole stress- 
energy tensor with the same eigenvalue e. r\ is the shear viscosity of the medium. 

The Navier-Stokes form is conceptually simple but introduces unphysical super- 
luminal signals that lead to numerical instabilities. 



The second-order Israel-Stewart formalism 36 - 38 avoids this super-luminal 



propagation, as do more recent approaches 39 40 . In the Israel-Stewart formal- 
ism for a conformal fluid 34 , the stress-energy tensor is decomposed as T^ v — 



Tff + tt^ . The evolution equations are 

9 M r^ = 0, (5) 

A£A> CT d (T 7r a ' 3 = -— (tt"" - S"") - \^{d a u a ) , (6) 
t ot 3 

where for brevity we neglected to show heat-flow, vorticity and terms that turn 
out to be numerically irrelevant. The role of vorticity in heavy-ion collisions when 
including fluctuations has also recently been studied [41] . 

We note that the Israel-Stewart formalism is based on two main choices: The 
14-moment approximation to truncate the single particle distribution function, and 
the choice of the equations of motion for the dissipative currents [40] , which is 
not unique. We will come back to this point when discussing recent theoretical 
developments in Section [6j 

Simulations of bulk dynamics in heavy-ion collisions using the Israel-Stewart 



formalism have been performed in 2+1 dimensions 14 31 42-45 , assuming boost 



invariance in the longitudinal direction. This is a reasonable assumption around 
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mid-rapidity for high enough energy collisions, particularly at LHC energies. 2+1 
dimensional simulations were also performed using the equivalent Ottinger-Grmela 
46 formalism 47 . More recently, 3+1 dimensional viscous calculations that are 



required away from mid-rapidity and at lower collision energies have also become 



available 48 -50 



Having established the standard theoretical framework, in the following we will 
discuss recent developments in the field of relativistic viscous hydrodynamics and 
its application to heavy-ion collisions. 



4. Equation of state 

Typically, the equation of state ([3]) used in hydrodynamic simulations of heavy- 
ion collisions is determined from lattice QCD calculations combined with a hadron 
gas model. A detailed comparison of different lattice equations of state has been 



performed 51 and it was found that the different equations of state 45 52-55 



which have different results for the trace anomaly and the speed of sound, lead 
to different evolution of the momentum anisotropy when used in an ideal hydro- 
dynamic simulation. However, the difference in the final spectra and elliptic flow 
results turns out to be negligible. This and uncertainties in the initial state make a 
precision determination of the equation of state from the comparison of hydrody- 
namic calculations to hadronic experimental data impossible. On the other hand, 
using the latest lattice equations of state within hydrodynamic calculations can po- 
tentially allow for the reliable extraction of medium properties such as transport 
coefficients from hadronic observables. It is particularly encouraging that lattice 
QCD calculations from the HotQCD and Wuppertal-Budapest collaborations now 
agree on the value of the (pseudo-)critical temperature for chiral symmetry restora- 



tion: T c = 154 + 8 (stat.) ±1 (sys.) MeV 56 vs. T c = 155 + 3 (stat.) ±3 (sys.) 



MeV 57 . Despite this agreement, systematic uncertainties remain in the interac- 
tion measure. Results from the HotQCD collaboration agree well with those from 
the Wuppertal-Budapest collaboration below T c , but are systematically larger by 
approximately 25% above T c . 

The production of thermal photons and dileptons is potentially more sensitive to 
the equation of state, because they are produced throughout the evolution and hence 



probe the dynamics of the system more directly 58 -60 . Another relevant aspect 



is the inclusion of chemical freeze-out in the equation of state (and the freeze-out 
procedure) to reproduce the correct final particle ratios [3| [6T|[62] . 

The determination of the equation of state from lattice QCD has thus far been 
restricted to zero baryon chemical potential, which is a good approximation at LHC 
and top RHIC energies at mid-rapidity. However, away from mid-rapidity and at 
smaller energies the finite net-baryon density becomes an important factor. Because 
lattice QCD is plagued by the sign-problem at non-zero baryon chemical potential, 
which prohibits the use of standard Monte Carlo techniques, the extraction of an 
equation of state at non-zero baryon chemical potential is very challenging. However, 
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at small values of the baryon chemical potential, the construction of a more general 
equation of state using a Taylor expansion of the pressure |63| is possible. Recently, 



such equation of state for finite net-baryon densities has been parametrized 64 



using the results of a lattice QCD analysis of the Taylor coefficients 65 66 



5. Bulk viscosity and limits of second order viscous hydrodynamics 

QCD is a non-conformal theory and thus contains a finite bulk viscosity. The effect 



of bulk viscosity on elliptic flow 67 and the combined effect of shear and bulk 
viscosity 68,69 have been studied. Bulk viscosity is expected to peak (possibly 
along with its relaxation time) around the critical temperature T c , where the system 
may develop large correlation lengths [9 70 - 72 . Within the range of applicability 



of second order viscous hydrodynamics, where viscous corrections to the thermal 
equilibrium distribution function remain small compared to the equilibrium part, 
corrections from bulk viscosity are found to be small compared to those from shear 
viscosity if the bulk relaxation time peaks around T c |69| . This means that if bulk 
and shear viscous corrections do not become so large as to render the hydrodynamic 
expansion invalid for the relevant part of a heavy-ion collision, the extraction of the 
shear viscosity should be possible to reasonable accuracy when neglecting the effect 
of bulk viscosity. 

Recent studies of the bulk viscosity have observed that while the bulk viscosity 
itself scales as the second power of conformality breaking ( ~ ri(c 2 s — 1/3) 2 , the 
correction to the distribution function Sf scales as the first power. [73] Corrections 
to the spectra are therefore dominated by viscous corrections to the distribution 
function, and reliable bounds on the bulk viscosity require accurate calculations 
of Sf in the hadronic phase. We will return to this point when discussing viscous 
corrections to the distribution function in Section |8.1| However, also in this study 
the effect of bulk viscosity on the p-r-integrated pion w 2 at RHIC energy is found to 
be small. 

Apart from viscous corrections to the thermal equilibrium distribution func- 
tions becoming large, bulk and shear viscous corrections can lead to a negative 
longitudinal pressure, argued to potentially cause the breakup of the system into 



droplets 74-76 . However, second order viscous hydrodynamics can only hint at 
when such cavitation could happen but is not suited to describe the process, as it 
happens outside the range of its applicability. 

6. Progress beyond second order viscous hydrodynamics 

The problem with the viscous correction to the stress energy tensor becoming larger 
than the equilibrium part, leading to negative pressure, arises especially in the early 
stage of the evolution, when the local momentum distribution is not yet equilibrated 
but highly anisotropic due to rapid longitudinal expansion. Recent progress has been 
made in describing this early time evolution and late time hydrodynamics within the 
same framework by performing the hydrodynamic expansion around an anisotropic 
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distribution 77-80 . This results in new equations of motion, including one for 
the degree of anisotropy of the distribution function. The procedure can reproduce 
both the limits of free streaming and ideal hydrodynamics and results in second 
order viscous hydrodynamics when expanding around a small anisotropy parameter, 
which so far has been shown in the one dimensional |77| and boost-invariant 2+1 
dimensional cases [8l] . Another interesting aspect of this anisotropic hydrodynamics 
is the fact that entropy production is modified as compared to second order viscous 
hydrodynamics. One would expect an increase in entropy production up to some 
value of viscosity over entropy density followed by a decrease as one approaches the 
limit of free streaming. This behavior is reproduced in anisotropic hydrodynamics 



77 but missed in second order viscous hydrodynamics as the free streaming limit 
is far outside its range of applicability. 

Another logical step is to expand to third order in gradients [82]. In the third 
order theory, numerical differences to second order Israel Stewart theory are com- 
pletely negligible for rj/s = 0.05 but become significant for rj/s > 0.2. It has also 
been pointed out that the hydrodynamic equations depend on the details of their 
derivation. While Israel and Stewart used the second moment of the Boltzmann 
equation to derive hydrodynamic equations for the dissipative currents, in recent 



work 83 the definition of the latter was used directly. This leads to equations of 
motion of the same form but with different coefficients. Microscopic transport cal- 
culations 84 , 85 show very good agreement with solutions to these equations of 



motion up to rj/s ~ 3, while the Israel-Stewart equations show large differences for 
rj/s > 0.2. Furthermore, solutions to the Israel-Stewart formalism for heat flow do 



not agree with direct numerical solutions of the Boltzmann equation 86 , even when 
the Knudsen number is very small (close to the ideal case). An improved deriva- 
tion of relativistic dissipative fluid-dynamics from kinetic theory without using the 



14-moment approximation was recently presented 87 88 and provides a good de- 
scription (compared to direct numerical solutions of the Boltzmann equation) of all 
dissipative phenomena, including heat flow. This demonstrates that the details of 
the derivation of the system of relativistic dissipative fluid-dynamic equations are 
relevant in particular when matching to kinetic theory at late times, as done in 



simulations using so called hadronic afterburners (see Section 8.2). 



7. Temperature dependent rj/s 



By comparison of experimental data with viscous hydrodynamic calculations using 
constant rj/s we can extract at best an effective (r)/s). In reality rj/s should depend 
on the local temperature of the medium, dropping from large values at high tem- 
peratures to a minimum around T c , and rising again with decreasing temperature 



in the hadronic phase F7j[89 ■ 92 



In fact, calculations with constant rj/s find good agreement with experimental 
data using a smaller effective value of rj/s at RHIC than at LHC 20 93 , in line 
with the expected increase of (rj/s)(T) at high temperatures. 
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It has been shown that when including a modeled temperature dependence in 
calculations at RHIC energies, the details of r]/s(T) in the quark-gluon plasma phase 
have little influence on the final elliptic flow result, while hadronic rj/ s{T) modifies 



V2 strongly 92 . At the highest LHC energies the conclusion is the opposite: weak 
dependence on the hadronic rj/s but strong dependence on rj/s{T) in the QGP 
phase. Interestingly, at RHIC energies, a significant dependence has been found on 
the minimum value of rj/s(T) around T c [94] . 

In addition, a strong dependence of the entropy production on the initial value 



of n 1 *" was found in 95 when starting with a large r)/s. After correcting for this to 
reproduce final particle spectra, differences in V2 were however minor. Yet, partic- 
ularly when including a temperature dependent f]/s, it is essential to gain a better 
understanding of the pre-equilibrium stage in heavy-ion collisions to determine the 
initial conditions for viscous hydrodynamics. 

Recent simulations have also incorporated a temperature dependent (j]/s)(T) 



20 , 50 and bulk viscosity over entropy density ratio (£/s)(T) [50] in 3+1 dimen- 



sional simulations. 



8. Conversion into particle spectra 

To compare to experimental measurements, the evolved stress-energy tensor has to 
be converted into particle spectra. It is possible to perform the kinetic freeze out 
in the hydrodynamic simulation, given a freeze-out condition. The typically used 
condition is the system reaching a certain freeze-out temperature or energy density. 
However, it has been pointed out 96 .97 that a realistic description should involve 
a comparison of the system's expansion rate with the particle scattering rate: when 
the scattering rate is not significantly larger than the expansion rate the system 
will freeze out. 



Particle spectra can be computed using the Cooper-Frye formula 98 



E™= JdZ M pPf(T,ppu' i ,*v), (7) 

where f^T^p^u^, tt^) is the particle distribution function, which is thermal in the 
ideal case and includes corrections proportional to 7r MI/ in the viscous case. E is the 
freeze-out surface, generally a three-dimensional surface in 3+1 dimensional space- 
time, which characterizes the distribution of time and space where the system freezes 
out. 

Another, perhaps more realistic way to deal with the late stage of a heavy- 
ion collisions when the system becomes dilute is to switch from a hydrodynamic 



description to a hadronic cascade simulation 99 - 105 . Conversion into individual 
particles can be done by sampling the spectra determined from ([T]) . These particles 
are then evolved in the hadronic cascade and freeze-out happens automatically as 
the system continues to expand and interactions become very rare. 

We discuss several important aspects of the conversion of hydrodynamic quan- 
tities into measurable particles in the remainder of this section. 



January 25, 2013 2:3 



HYDRODYNAMIC MODELING OF HEAVY-ION COLLISIONS 11 



8.1. Viscous corrections to particle distribution functions 

As alluded to above, when translating the dissipative T^ v to particles in the Cooper- 
Frye formalism (98 , corrections to the distribution function Sf have to be taken into 
account: 

ii&» = 5>/ |^9f (/on+<%) (8) 

for an N component system with d n the degeneracy of species n. The ansatz 

Sfn = §kfo(l±h)p a fx(p)^f, (9) 

where p a is a unit vector in the a direction, leaves some freedom and the usual 
procedure is to assume that all coefficients C n are equal, even though they should 
depend on the individual particle species' interaction rate 106 , and use xip) = P 2 > 



which is derived within a relaxation time plus Boltzmann approximation: 

Sf n = / 0n (l ± /on)pV^ 2 r^ 2 V n . (10) 

It has been shown 107 that x(p) °c P a ^ where a can take on values from 1 to 2, 
which is the case for example for a system with radiative and elastic energy loss 
that has x(p) °c p 1 ' 38 - The exact form of the correction has a sizable effect on the 
Pt differential elliptic flow for p^ > 1 GeV. It should be noted, however, that the 



analysis of experimental data 108 indicates x(p) ot p 2 . A more general problem 
is that corrections Sf can become large compared to /o- For hadrons this problem 
can be reduced when using a hadronic afterburner and switching at intermediate 
temperatures of ~ 160 Me V, but for photons that are produced throughout the 



whole evolution this becomes a serious concern 109 . 

As mentioned earlier, corrections from bulk viscosity are most prominent in 
modifications of 5f, which is proportional to the conformal breaking. Therefore 
accurate calculations of Sf in the hadron resonance gas are required. 73 While 
effects of bulk viscosity on integrated seem small, bulk viscous effects need more 
detailed studies to determine £/ ' s and rj/s as functions of temperature in the future. 



8.2. Hadronic afterburner 

Hydrodynamic description of a system applies when the mean free paths of con- 
stituent particles are small compared to any macroscopic length scale. It is then well 
justified to use hydrodynamics for the initial phase of the QGP evolution since the 
density is so high. As the system expands and cools, the mean free paths decrease 
and the ratio r//s increases. Eventually, when the density becomes low enough, the 
fact that each species has a different mean free path begins to matter. This is be- 
cause differing mean free path not only implies different speed of response to the 
change in the collective flow, but it also implies that the collective description is 
becoming inadequate. On the other hand, the kinetic theory model of heavy ion 
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collisions works best when the system is dilute and the mean free paths long. Fur- 
thermore, it is only natural that the mean free path of each species is different in 
the kinetic theory models. 

If one can assume that there is a range of temperature and density that the 
validity of both the hydrodynamic approach and the kinetic theory approach over- 
lap, then it becomes natural that the hydrodynamic description should give way to 
a kinetic theory description when temperature and density fall within this range. 
Although it is not rigorously proven that such a transition region exists, evidence 
shows that including such transition in the simulation does improve the agreement 
between theory and experiments 33,93,99, 100, 102 |l 10 - 113 , especially in the bary- 



onic sector. 

Coupling between an event-by-event 3+1D viscous hydrodynamics and a ki- 



netic theory model has been recently accomplished 114 . The publicly available 
UrQMD code (v. 3.3pl) [TlO] was joined with the event-by-event 3+ ID viscous hy- 
drodynamic simulation MUSIC, taking special care in matching the rj — r coordinate 
system used in MUSIC to the t — z system used in UrQMD. To couple MUSIC to 
UrQMD, first the isothermal 3-D hyper-surface at T = 170 McV is determined in 
the hydrodynamic simulation. Next the transition hyper-surface is sampled using 
the Cooper-Frye formula Eq. Q to populate hadrons in i — z space. These hadrons 
are then propagated backward in time without collisions to the common UrQMD 
initial time to- 

At this point, the regular UrQMD evolution takes over, but with the important 
difference that the particles are allowed to interact only after they emerge from 
the isothermal hyper-surface. This is done because different cells reach the tran- 
sition temperature at different times whereas UrQMD needs to evolve particles in 
Minkowski time. In this way, particles are simulated as emerging from the plasma 
at different times. 



Three main results from this study are shown in Figs. [2pj and g] Hybrid cal 



culations with microscopic cascade simulations typically use oversampling to save 
computing time. In the case of MUSIC with UrQMD in each centrality bin, 100 MU- 
SIC events are sampled 100 times each for RHIC and 1000 MUSIC events are sampled 
10 times each for LHC. Light hadron spectra for 7r~, K~ and p for RHIC are shown 
in the left panel of Fig. [2j Experimental data from PHENIX 115 is reasonably well 



reproduced for the shown — 5 % centrality bin. For all three species, the descrip- 
tion of the data is either improved over or of about the same quality as the pure 
hydrodynamics results. In the second panel, LHC predictions are compared with 



the ALICE data shown at QM2012 116 



For V2(pt) shown in Fig. [3j the p result is much improved with the afterburner 
at RHIC. This may be due to the better description of the finite baryon mean free 
path in UrQMD compared to hydrodynamics. For the LHC, one may argue that the 
afterburner improves the description a little, but the statistical error at this point 
is too big to make a definite statement. Calculations with afterburner lead to larger 
v z{pt) m general (see Fig. [4]). However, this needs more careful study to quantify. 
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Fig. 2. Transverse momentum distribution of identified particles at RHIC and LHC for 0-5% 
most central collisions compared to MUSlC+UrQMD calculations. RHIC result is based on 10,000 
events (100 UrQMD events on each of 100 hydro events). LHC result is based on 1,000 events (10 
UrQMD events on each of 100 hydro events). For RHIC, the results of pure MUSIC events are also 
shown. Experimental data from PHENIX |1 15] and ALICE (preliminary) [116| . 




Fig. 3. The dependent V2 from RHIC and the LHC in the 10 — 20 % centrality class compared 
to MUSIC and MUSlC+UrQMD calculations. RHIC result is based on 10,000 events and LHC result 
is based on 1,000 events. Experimental data from STAR [117] and ATLAS [118] , 
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Fig. 4. The pj> dependent triangular flow from RHIC and the LHC in the 10 — 20 % centrality 
class compared to the MUSIC and MUSlC+UrQMD calculations. RHIC result is based on 10,000 
events and LHC result is based on 1,000 events. Experimental data from PHENIX [119| and 
ALICE [120]. 
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9. Event-by-event hydrodynamics 
9.1. Overview 

It is now generally accepted that fluctuating initial conditions for hydrodynamic 
simulations of heavy- ion collisions are essential for the exact determination of collec- 
tive flow observables and to describe features of multi-particle correlation measure- 
ments in heavy-ion collisions 17 18 48 49 121 ■ 136 . Real event-by-event hydrody- 



namic simulations have been performed and show modifications to spectra and flow 
from "single-shot" hydrodynamics with averaged initial conditions 48 , 49 , 134 136 



An important advantage of event-by-event hydrodynamic calculations is the possi- 
bility to study higher flow harmonics such as ^3, which is entirely due to fluctuations 
like all odd harmonics. Different v n depend differently on rj/s and the details of the 
initial condition, which is determined by the dynamics and fluctuations of partons 
in the incoming nuclear wave functions. This observation can be used to deter- 
mine these long sought after details of the initial state and medium properties in 
heavy- ion collisions by performing a systematic analysis of all harmonics v ni up 
to e.g. n — 6 as a function of 77/s and the initial state properties and compare 
to experimental data. First predictions of ^3 from hydrodynamic simulations 48 
agree extremely well with experimental data from RHIC |137| . Furthermore, it has 
been shown that at low p T (and |Ar?| > 1 (ALICE), \Arj\ > 2 (ATLAS)), the main 
features of dihadron correlations in the angular difference A0 between the hadron 
momenta can be described by flow, i.e., the sum of v\ to vq only |138 139| . The 
double-peak structure on the away-side is hence described mostly by (triangular) 
flow as some works had predicted 17 131 . Detailed measurements and calculations 



of new quantities, such as event plane correlations or event-by-event distributions 
of flow coefficients have appeared over the last year and will be reviewed in this 
section. 



9.2. Initial state models 

There are several ways for describing and generating fluctuating initial conditions 
for hydrodynamics. Generally, this means the energy density distribution at the 
"thermalization time" and in some cases the distribution of flow velocities. One of 



the most commonly used models is the Monte-Carlo (MC) Glauber model 140 
In its simplest implementation uncorrelated nucleons are sampled from measured 
density distributions. Then the two nuclei are arranged according to a random im- 
pact parameter b and projected onto the transverse x-y plane, assuming straight 
line trajectories for all nucleons. Interaction probabilities are then computed using 
the relative distance between two nuclei and the measured nucleon-nucleon inelastic 
cross section. Every wounded nucleon (sometimes also a fraction of binary collision 
points) is then assigned an energy or entropy density, parametrized e.g. as a two 
dimensional Gaussian in the transverse plane. The model can be improved by the 
inclusion of many-body correlations between the nucleons. However, it was shown 
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that the inclusion of realistically correlated configurations does not modify the ge- 
ometry, characterized by the dipole asymmetry and triangularity of the distribution, 
compared to the uncorrelated case |141| . 

The MC-KLN (Monte-Carlo Kharzeev-Levin-Nardi) model, another commonly 
used initial state model, is based on the color-glass-condensate framework. Within 
this framework, one takes into account the well established feature of QCD that at 
small Bjorken-x, a novel regime governed by large gluon densities and non-linear 
coherence phenomena takes over 142 . These high gluon densities correspond to 



strong classical fields, permitting calculations of the wave function using classical 
techniques. Quantum corrections are then incorporated via non-linear renormal- 
ization group equations such as the JIMWLK or in the large N c limit the BK 
equations 143 144 that describe the evolution of the wavefunction towards smaller 



x. Certain saturation models incorporate this evolution via parametrizations of the 
unintegrated gluon distribution or other quantities and we describe some of these 
approaches in detail below. 

The MC-KLN model in particular is built upon the KLN model |145p46] , which 
was improved to the so called f(actorized)KLN model and extended to a full Montc- 
147p48 



Carlo version 



The main ingredient for this description is the fcy-factorized 
expression for the gluon multiplicity, essentially a convolution of the unintegrated 
gluon distributions 4>a/b of the two colliding nuclei: 



(Prxdy 



4iV e 
N2- 



Pt 



d 2 p T f PT d 2 k T 



Pt 

x <f> A ( Xl , (p T + k T ) 2 /4)0 B (x 2 , (p T - k T ) 2 /4) (11) 



Here, px and y denote the transverse momentum and the rapidity of the pro- 
duced gluons, respectively, r-r is the position in the transverse plane. The light- 
cone momentum fractions of the colliding gluon ladders are then given by x\^% = 
Pt exp(±y)/y / s, where s denotes the center of mass energy. 

A few remarks on fc^-factorization are in order. A general shortcoming of kr- 
factorization is that multiple scatterings for px < Q s are omitted even though they 
should be present because occupation numbers are high. So the approach is strictly 



not valid in the fully non-linear case of A+A collisions 149 . In fact, Eq.( 11 ) cannot 



be derived for collisions of two dense systems but is rather an extrapolation from 
the result in the "dilute" p+A limit. The cutoff of the kr integral at pr, used in 
the MC-KLN model, is introduced by hand to make the spectrum infrared finite. 
Different implementations of fcx-factorization and comparisons to more appropriate 



classical Yang-Mills calculations have been discussed in the literature 150 



The KLN approach uses a particular parametriza tion of th e unintegrated gluon 
distribution, depending on the saturation scale Q^ B 145|l46 . Apart from the men- 
tioned problems with ^-factorization for A+A collisions, there are several problems 
with the original KLN approach, in particular that the saturation scale in one nu- 



cleus depends on properties of both nuclei 151 , and that the limit at the edge of a 



nucleus does not approach the unintegrated gluon distribution of a single nucleon. 
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These problems are cured in the fKLN approach |147| where the result for the gluon 
multiplicity is truly factorized (unintegrated gluon distributions of nucleus A (B) 
only depend on properties of nucleus A (B)) 



dN g f d 2 p T 



cPrxdy J p\, 



d 2 kTPAPB(t>A(T A /p A )4>B{T B /pB) , (12) 



where pa/b is the probability for finding at least one nucleon at a given transverse 
coordinate in nucleus A/B, with T^/b being the thickness function. The implemen- 
tation of Monte-Carlo sampling also takes care of correct behavior at the edges. The 
initial energy density is obtained by computing the transverse energy distribution 
of the produced gluons, i.e., by including another factor of px in the integrals in 
Eqs. ( 11|12 ). The model has been further extended to include running-coupling BK 



(Balitsky-Kovchegov) evolution of the gluon distributions in the rcBK model 152 
MC-KLN-models, as the MC-Glauber model, lack negative binomial fluctuations 
of multiplicities for a given number of participants. In a recent extension, these 



fluctuations have been added to the MC-KLN model by hand 153 . In this imple- 
mentation the correct correlations in the transverse plane are obtained by choosing 
the grid cell size to be 1/Q S . 

A more recent improved color-glass-condensate based approach is the IP-Glasma 
model [154 155] , which combines the IP-Sat (Impact Parameter dependent Satu- 



ration Model) model 156 157] of high energy nucleon (and nuclear) wavefunctions 



with the classical Yang-Mills (CYM) dynamics of the glasma fields produced in a 
heavy- ion collision 158 - 163] . After fixing the free parameters of the IP-Sat model 



by fits to small x HERA deeply inelastic scattering (DIS) data off protons and fixed 



target nuclear DIS data 164 165 , the IP-Sat model provides an excellent descrip- 
tion of these data. The IP-Glasma model includes fluctuations of nucleon positions 
as well as sub-nucleonic fluctuations of color charges, a feature missing in most 
other initial state models. Another advantage is that the model does not rely on 
fey-factorization, which is strictly only valid when at least one of the sources is dilute 
(as in p + p and p + A collisions) as discussed above. Furthermore, the IP-Glasma 
includes the non-linear pre-equilibrium evolution of the initial gluon fields. This 
leads to the build-up of initial flow and an independence of the exact time when 



one switches to hydrodynamics 20 . The early stage dynamics is however not fully 
included. Instabilities triggered by quantum fluctuations, and subsequent strong 
scattering of over-occupied fields, may lead to rapid isotropization and quenching 
of IP" to reasonable values justifying the use of viscous hydrodynamics already 
at early times. These unstable dynamics require a full 3+1 dimensional simulation 
including a realistic description of quantum fluctuations, which has not yet been 



fully achieved. However, significant progress is being made 166 ■ 170 and the IP- 
Glasma model can be extended to include these important effects. Matching of the 
full stress-energy tensor, including viscous corrections and flow, to the hydrody- 
namic simulation will then be possible. Another possibility is to couple this initial 
condition to anisotropic hydrodynamic simulations 77 - 80 described in Section [6] 
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The additional color charge fluctuations in the IP-Glasma model naturally lead to 
negative binomial fluctuations in the event-by-event multiplicity and the correla- 
tion length of the fluctuations in the transverse plane is of the order of the inverse 
saturation scale 1/Q S as desired. 

In Fig. [| we show a comparison of initial energy densities from an MC-Glauber, 
the MC-KLN and the IP-Glasma model using the same distribution of nucleons 
in the incoming nuclei. In the MC-Glauber model every wounded nucleon was as- 
signed a two dimensional Gaussian energy density with a width of <7o = 0.4 fm. The 
MC-KLN result was obtained using the publicly available code mckln-3.52 |171| . 
IP-Glasma results are shown for two different times, r = 0.01 fm/c and r = 0.2fm/c 
after Yang-Mills evolution. The evolution smoothens the initially very distinct struc- 
tures noticeably. Because of the additional subnucleonic fluctuations, the IP-Glasma 
model produces the finest granularity, typically leading to larger fluctuation driven 



odd eccentricities 154 155 





Fig. 5. Comparison of the initial energy density (arbitrary units) produced by the MC-Glauber, 
MC-KLN and IP-Glasma models. All events have the same configuration of nucleons and impact 
parameter 6 = 4 fm to emphasize how different model descriptions affect the structure of the energy 
density. The finest structure is obtained in the IP-Glasma model, which includes subnucleonic color 
charge fluctuations. Yang-Mills evolution to r = 0.2 fm/c smoothens this structure before it enters 
a hydrodynamic simulation. 



Apart from MC-Glauber and CGC based frameworks, there are several parton- 
and hadron-cascade models that are being used to determine fluctuating initial 
conditions. These are for example UrQMD [135], EPOS [172], and AMPT [T73"1|T74 



all using Monte-Carlo techniques to compute initial particle production and then 
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converting the soft part of the spectrum into the bulk energy density distribution 
used in hydrodynamic simulations. They also provide initial flow and in principle 
the full stress-energy tensor including viscous corrections. 



9.3. Hydrodynamic fluctuations 

In addition to fluctuations in the initial state, there should be fluctuations occurring 
during the hydrodynamic evolution, simply due to the finite number of particles 
present in a real heavy-ion collision. These hydrodynamic fluctuations 175||176 are 
directly related to the hydrodynamic properties of the matter, in particular, owing to 
the fluctuation-dissipation theorem, their amplitude is governed by the viscosities. 
This offers the possibility to constrain the viscosity of the strongly coupled quark- 
gluon plasma independently from the traditional analysis of anisotropic flow. Full 
simulations including this type of fluctuations are still to be developed. 



9.4. Numerical schemes 



When dealing with fluctuations and trying to learn about the physical viscosity 
of the system, one needs numerical schemes to solve the hydrodynamic equations 
that can deal well with large gradients and have a small artificial viscosity. It has 
been shown that in the case of heavy-ion collisions, several algorithms do a good 



job in this respect 177, 178 . These are, in particular, the SHASTA algorithm 179 



which after finding a low-order solution with large numerical diffusion, corrects the 
result using anti-diffusion fluxes in a second step, and the Kurganov-Tadmor (KT) 



algorithm 180 , an improvement of the Nessyahu-Tadmor algorithm (NT) [181], 



which is Ricmann-solver free and can be viewed as an extension of the Lax-Friedrichs 
(LxF) scheme. Recently an algorithm with a relativistic Riemann solver geared 



towards simulating heavy-ion systems 177 has also been developed. Comparisons 



of the performance of the first three types of algorithms in solving the relativistic 



Riemann problem were performed 178 and no clear difference was found in either 



performance or accuracy of the algorithms at a given lattice size. The analysis was 



extended 177 to include the Riemann solver based algorithm, which naturally is 



found to have some advantage in solving the Riemann problem in the test scenario. 



9.5. Recent results and comparison to experimental data 

Recently, event-by-event hydrodynamic calculations have produced a plethora of 
predictions and explanations of a wide range of experimental data, for which the 
inclusion of initial state fluctuations is the essential ingredient. Not only do these 
simulations allow for the determination of all average higher flow harmonics beyond 
elliptic flow, but further details like the event-by-event probability distributions of 
the flow harmonics or event plane correlations can be computed and compared to 
experimental data. 

Simulations using the IP-Glasma initial state model have been particularly suc- 
cessful in describing both the px dependent and integrated v n at both RHIC and 
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LHC energies |20| . The agreement with experimental results from LHC shown in 
Fig. [6] is particularly striking. 
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Fig. 6. Left: Root-mean-square anisotropic flow coefficients (i) 2 } 1 / 2 in the IP-Glasma model |2p| , 
computed as a function of centrality, compared to experimental data of v n {2}, n £ {2,3,4}, 
by the ALICE collaboration [182] (points). Right: Root-mean-square anisotropic flow coefficients 
(dJ) 1 / 2 as a function of transverse momentum, compared to experimental data by the ATLAS 
collaboration using the event plane (EP) method [22] (points). Bands indicate statistical errors. 



This agreement indicates that initial state fluctuations in the deposited energy 
density, translated by hydrodynamic evolution into anisotropies in the particle pro- 
duction, are the main ingredient to explain the measured flow coefficients. 

Because of this feature, some effort has been concentrated on characterizing the 
initial state in a way that ties it directly to the measured flow. The simplest way of 
doing so is to compare the initial eccentricities of the system 



■\/{r n cos {ruj))) 2 + (r n sin(n</>)) 2 



(13) 



to the final flow harmonics v n . However, in particular for and higher harmonics, 
the nonlinear nature of hydrodynamics becomes important |183| and more accurate 
predictors for flow coefficients involve both linear and nonlinear terms, e.g. v§ has 



contributions from £5 and £263, and it was shown 184 that the nonlinear term 
becomes more dominant with increasing viscosity. 

The fact that linear terms are damped more by viscosity leads to a growing 
correlation of different event planes 



1 (sin(n</>)) 

arctan —. -rr , 
n (cos(n<p)) 



(14) 



with increasing viscosity 184 , a result that is in line with findings in a different 



work [185], where experimental data on event plane correlations from the ATLAS 



collaboration 186 was compared to hydrodynamic calculations in different scenar- 
ios. 
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It has been argued that the most precise value for the shear viscosity to entropy 
density ratio r)/s of the quark-gluon plasma can be obtained by a fit to 2?T-hitegrated 
v n measurements in ultra-central collisions, because in those fluctuations dominate 
all v n and uncertainties in £2, which varies most from one initial state model to the 



next, are minimal. 10 At the very least, the study of ultra-central collisions can 
constrain the possible fluctuations implemented in different models by comparison 
to experimental data. 

Another way to constrain initial state models with fluctuations is the study of 
event-by-event distributions of flow harmonics v n , which have recently be deter- 



mined experimentally [211. It was found 187 that these distributions are almost 



independent of the details of the hydrodynamic evolution, like the shear viscosity 
to entropy density ratio. In fact, the distributions of initial eccentricities provide 
already an excellent approximation of the measured v n distributions when scaled by 
the mean value. This allows for an almost independent determination of the initial 
state model and the transport properties of the evolving system. As shown in Fig. 
[7] the IP-Glasma model provides a very good description of these probability distri- 



butions 20 and differences between the e„ and v n distributions are only visible in 
the tail where nonlinear effects in the hydrodynamic evolution are potentially more 
important. 

While most simulations take into account fluctuations of the energy density 
and flow velocities in the transverse plane only, some studies have also considered 



fluctuations in the longitudinal direction 41 174 189], finding a reduction of pion 



elliptic flow compared to calculations with only fluctuations in the transverse plane 



174 . This is an interesting possibility and needs to be investigated more closely in 
the future. 

Event-by-event hydrodynamic calculations have also been utilized to compute 
dihadron correlations in heavy-ion collisions. When including local charge conser- 
vation, agreement with the two-dimensional two-particle correlation data in rela- 
tive azimuthal angle and pseudorapidity at soft transverse momenta (px < 2 GeV) 



is found 190 . It should however be noted that the origin of the almost boost- 
invariant initial state, necessary to reproduce the measured "ridge" structure in 
pseudo-rapidity is not described by hydrodynamics - it is merely maintained by the 
hydrodynamic evolution, and requires an independent explanation. 

Recent efforts have pushed these studies further to include p + A and even 



p + p collisions 191 ■ 193 . Whether the system size and viscous corrections in such 
small systems allow for the application of hydrodynamics is however questionable. 
Saturation models have reproduced correlation data over a wide kinematic range, 



194- 196 showing that additional hydrodynamic flow is in fact required to describe 



the data from heavy-ion collisions, but has no room in smaller systems like p+p 194 . 
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Fig. 7. Scaled distributions of 1)2, f3 and 114 as well as £2, £3 and £4 from the IP-Glasma model [20| 
compared to experimental data from the ATLAS collaboration |21|188] . Using 750 (0-5%) and 1300 
(20-25%) events. Bands are systematic experimental errors. 



10. Summary, Conclusions and Outlook 

Relativistic viscous hydrodynamics has been extremely successful in describing the 
bulk properties of heavy-ion collisions at RHIC and LHC, ranging from particle 
spectra to anisotropic flow and correlations. In particular the use of event-by-event 
hydrodynamic simulations together with models for the fluctuating initial state 
has dramatically increased the amount of successful predictions, such as the dis- 
tributions of higher harmonics v n , correlations of event planes, etc. Comparison of 
experimental data to hydrodynamic simulations thus allows to extract properties 
of the matter created in heavy-ion collisions, such as the fact that the system is 
strongly interacting, and more quantitative measures like the shear viscosity to en- 
tropy density ratio rj/s. Hydrodynamics is at this point in time the best tool for the 
determination of such fundamental properties of a hot and dense quantum-chromo 
dynamic system. 

Despite these major successes, several unknowns and uncertainties remain. In 
particular it is not understood how the system reaches a thcrmalized state that is 
necessary for the usual hydrodynamic framework to be applicable. In fact, it is not 
settled whether the system actually does thermalize by the time we start applying 
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hydrodynamics or whether momentum distributions remain largely anisotropic up 
to late times and hydrodynamics should be modified to accommodate this situation. 

Related to this is the question of how to initialize viscous corrections and whether 
they become too large in fluctuating events with large gradients. This would take us 
out of the domain of validity of hydrodynamics. The issue is particularly problematic 
when studying temperature dependent rj/s that can become very large at large 
(T ^> T c ) and small (T < T c ) temperatures. Temperature dependent values of 
additional transport coefficients such as bulk viscosity and relaxation times pose 
another large uncertainty that demands further detailed studies. 

In addition there remain uncertainties in the freeze-out description, be it a simple 
Cooper-Frye description with resonance decays or coupling to a hadronic rescatter- 
ing simulation. These lie for example in the determination and implementation of 
the freeze-out condition and the prescription of the conversion of energy densities 
to particle degrees of freedom. 

Current research is addressing these and other caveats and will hopefully shed 
further light on the complex details of the strongly interacting system created in 
heavy-ion collisions, and hence the properties of quant um-chromo dynamics under 
extreme conditions. 
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